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Abstract: Using the Euler-Maclaurin summation we calculate analytically the internal en- 
ergy for non-interacting bosons confined within a harmonic oscillator potential. The specific heat 
shows a sharp A-like peak indicating a condensation into the ground state at a well-defined tran- 
sition temperature. Full agreement is obtained with direct numerical calculation of the same 
quantities. When the number of trapped particles is very large and at temperatures near and 
above the transition temperature, the results also agree with previous approximate calculations. 
At extremely low temperatures both the specific heat and the number of particles excited from the 
condensate are exponentially suppressed. 
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Bose-Einstein condensation has now been experimentally demonstrated in magnetic 
traps of rubidium lithium Q and most recently sodium ||] gases. To a good approxima- 
tion one describes the trapping potentials as 3-dimensional anisotropic oscillators which 
in the rubidium experiments ||] have frequencies typically around u = 500 — 1000 Hz. For 
the sake of simplification we will here ignore the anisotropy and the weak interactions 
between the alkali gas atoms. The energy levels of one particle are then simply given as 
e n = Tiujn where the quantum number n takes the values n = 0, 1, 2 . . . when we drop the 
zero-point energy. Since each energy level has a degeneracy g n = (n + l)(n + 2)/2, the 
total number N of particles in such a trap at temperature T and chemical potential [i is 
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given by the Bose-Einstein distribution as 

9n 



N= y!iiL — (i) 

n=0 

where j3 = 1/ksT. The ground state has quantum number n = and thus contains 
iVo = A/(l — A) particles where A = exp (/?//) is the fugacity. We then have 

N = T^x +Ne (2) 

where the number of particles in the higher states is 

°° 9nXe- bn 



1 ~ ^ 1 - Xe~ bn (3) 



n=l 
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It depends only on the effective fugacity A = Aexp (—6) where b = hiv/k^T '. 



When the temperature is lowered, the fugacity A — > 1 and a finite fraction of particles 
starts to condense into the ground state. Since the total number of particles N 1, 
this will happen when the variable 6^1. In this region and at higher temperatures one 
can then approximate the sum in (||) by an integral. Including only the leading term 
~ n 2 in the degeneracy, one then has the semiclassical limit || . Defining the transition 
temperature to be the temperature at which the fugacity takes the value A = 1, it is then 
found |Q to be given by 

The resulting specific heat exhibits a sharp jump with the shape of a A. This is in contrast 
with the gas in zero potential where the specific heat is continuous. As a check, the 
partition function was also calculated using numerical summation. 

A more accurate approximation of the sum has been given by Grossmann and Holthaus 
0. Based upon the two leading terms in the degeneracy, they constructed a continuous 
density of states which then made it possible to approximate the sum by two integrals. 
This caused a small shift in the transition temperature by a factor ~ iV -1 / 3 (||) and a 
corresponding small rounding-off of the sharp peak in the specific heat when the particle 
number is not too large. In the experiments under consideration, N takes typically 
values in the range from 10 3 to 10 6 . 

Recently, these approximation methods in which the sum is replaced by integrals, have 
been criticized by Kirsten and Toms0. Instead, they propose to evaluate the sum directly 
by contour integration. In turns out, however, that this method is difficult to use in the 
transition region where they are forced to fall back upon a direct, numerical summation. 
Their results are thus quite similar to those of Grossmann and Holthaus ||. 

The standard method for summing a series like (H|) is Euler-Maclaurin summation ||. 
One then has 

n=b fX= b -. -, 

£/(n)= /_ dxf(x) + -[f(b) + f(a)} + -[f(b) -/»] + ••• (5) 

n=a Jx—a 

In our case the contributions from the upper limit will vanish. When T 3> hto/kB we can 
also safely ignore the contributions coming from the higher order derivatives at the lower 
integration limit. The sum (||) can then be written as 

N e = G 3 (A) + \g 2 {\) + G t (X) + ^j^j (| + j^-Sj (6) 

when we only include terms up to the first derivative. We have here introduced the new 
functions 

^ / T \ 1 f°° j xP\e~ bx 
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In the case of Bose-Einstein condensation in zero potential, the lower limit would have 
been x = and the functions would be equal to the g p (X) functions^]. These are equal to 
the more standard polylogarithmic functions 
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n=l 



We see that Li\{z) = — ln(l — z). Furthermore we will need Li p (l) = C(p)- By a partial 
integration one finds that the functions ([?]) can be obtained from the simple recursion 
relation 

G P+1 (X) = ^LhW + l£ yG p (A) (9) 

Since derivatives of polylogarithmic functions satisfy Lz p ' +1 (z) = Li p (z)/z, the new func- 
tions (0) can all be expressed by these. With the result (^) for the number of excited 
particles we then have 

N = T ^j + ^Li 3 (Xe- b ) + ^Li 2 (Xe- b ) + ^Li 1 (Xe- b ) 



1 Xe- b /31 

+ 4l3A^U + l3A^) (10) 



This equation determines the critical behaviour of the system. 

At temperatures below the semiclassical transition temperature (|j) most of the par- 
ticles are in the condensate consisting of Nq particles in the ground state of the harmonic 
oscillator. The fugacity A ~ 1 follows from 

A = ^L_ (11) 
N + l V ; 

when Nq S> 1. For higher temperatures we must use the full equation ( |l0|) to calculate A. 
The result is shown in Fig. | for N = 2000 and N = 20 000 particles. We see that when 
the number of particles gets to be very large, the transition into the ground state marked 
by having A = 1, becomes correspondingly sharp. In contrast, the effective fugacity A < 1 
at all temperatures as is clearly seen from Fig. [2[ The number of excited particles and 
as we will see, also the specific heat will thus be exponentially suppressed at very low 
temperatures. 

Above the transition temperature (|j) the two fugacities can be taken to be equal to a 
very good accuracy. This corresponds to having the parameter b <C 1 in this temperature 
region. The leading terms in the number of excited particles N e are the second and third 



terms in (10), i.e. 



N e = ±Li 3 (Xe~ b ) + A Li2 ( Ae ~ b ) + Q(l/b) (12) 
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Figure 1: The fugacity A as a function of T/T for N = 2000 and TV = 20000. A is close to unity at low 
temperature and starts to fall off rapidly around T — Tq. As TV is increased, the transition sharpens and 
the curve approaches the semiclassical result. 



Expanding the polylogarithms around b = we then get 

N e = pLi 3 (A) + ^Li 2 (\) + 0(1/6) (13) 

The two first terms here are the same as obtained in the calculation of Grossmann and 
Holthausg. They were obtained from a density of states based upon only the first two 
terms in the quantum mechanical degeneracy g n = (n 2 + 3n + l)/2. Their results are in 
surprisingly good agreement with the exact numerical results near and above the transition 
temperature. This is to a large extent explained by an apparent compensation of their 
neglect of the third term in the degeneracy by taking the integration from n = instead 
of from n = 1. Their results thus depend on the ordinary fugacity A instead of the 
effective fugacity A which appears in our more accurate analysis. As we have seen, it is 
the exponentially damped effective fugacity which determines the thermodynamics at very 
low temperatures. 

In the transition region we can use the approximate result (|l^) for the number of ex- 
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Figure 2: The effective fugacity A as a function of T/T for N = 2000 and N = 20000. A is always less 
than unity and vanishes at absolute zero. The maximum reached near the transition is seen to sharpen as 
N is increased. 



cited particles. The number of particles in the condensate thus varies with the temperature 
as 



^ - ( kBT \ 3 



Nn = N 



C(3) 



(14) 



Defining the transition temperature T c where Nq = we see that is is given by 



3C(2) (k B T c V 
2N \ Huj J 

\ C(2) 1 

2^2/3^1/3 _ 



1/3 



as obtained by Grossmann and HolthausQ. When the number of particles in the trap 
becomes very large, the transition temperature approaches the semiclassical result (|j). 
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Similarly, we get for the variation of the condensate ( |i~4| ) near the transition temperature 
No , fT\ 3 3C(2) 1 /T x2 



N 



1 - 1 5b 



2C(3) 2/3 VTq 



(15) 



This result is only approximately correct near the transition region where the effective 
fugacity A ~ 1. A more accurate result is obtained by numerically solving ( |To|) for A. 



0.18 



0.17 



0.16 



Exact sum 
Li 3 + Li 2 
Li 3 + Li 2 + Lij 
+ 5th term 
+ 6th term 



0.15 — 

0.910 0.915 0.920 

T/T 

Figure 3: The condensate ratio No/N as a function of T/T for N = 20000 just below the transition. 
The approximate results are seen to approach the exact numerical results as more terms in Eq. ( p"(i| ) are 
included. When the third derivative in the Euler-Maclaurin formula is included, a numerical agreement to 
five decimal places is obtained. 

The result is shown in Fig. || which gives the variation of the condensate just below 
the transition temperature. We see there the limited accuracy of just including the first 
two terms in N e when comparing to the exact result obtained by numerical summation. 
Using the Euler-Maclaurin formula up to and including the third-derivative term, we find 
a numerical agreement to five decimal digits. In Fig. || we show the condensate varying 
with the temperature down to absolute zero. As T — > the exponential suppression of 
excited particles becomes stronger and stronger. In fact, when T < huj/k, essentially all 
the excited particles will be in the n = 1 energy level of the oscillator. We then have 

N = N - 3e- huj/kBT (16) 
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when we neglect more exponentially suppressed terms. 
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Figure 4: The condensate ratio No/N as a function of T '/To for N = 2000 and N = 20000. iVo 
approaches N exponentially at extremely low temperatures. Near T = To the sharp transition found 
with the semiclassical approximation is rounded off by higher order corrections. In the figure one can not 
distinguish between the analytical and numerical results. 



The specific heat is obtained from the internal energy 



n=0 



,J3(E n -p) _ 1 



(17) 



In the condensate each particle has zero energy so that the sum really starts at the first 
excited level with n = 1. Again using the Euler-Maclaurin formula (|5|) we obtain 



U_ 



3 Li 4 (Xe- b ) + ^Li 3 (\e- b ) 



+ 



6 4 

1 \e~ b 



+ ^Li 2 (Xe- b ) + ^Li 1 (Xe- b ) 



4 1 - Ae" 



25 



Ae- 



(18) 



The specific heat can now be obtained as a function of the fugacity by taking the partial 
derivative with respect to temperature. It will depend on the partial derivative (dX/dT) 
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which can be obtained from (1C). The result is shown in Fig. [5] for different numbers 
of particles in the trap together with exact numerical results. At extremely low temper- 
atures the number of excited particles is exponentially small which will thus result in a 
correspondingly suppressed specific heat. 




Figure 5: The specific heat Cv/N as a function of T/T for N = 2000 and N = 20000. The discontinuity 
in the specific heat found with the semiclassical approximation is smoothened by the corrections. The 
sharpness of the peak is seen to increase with TV. 



When the number of trapped particles is very large, the first term in (18) gives the 
dominant contribution to the result. As a rough approximation, we can keep only this 
and setting A ~ 1 just below the transition temperature, we have 

i 



hu \ Tiuj J 



C(4) (19) 



The specific heat at the transition temperature is thus approximately 



CV (20) 



Nk B C(3) 

i.e. almost a factor four larger than for a classical gas in the oscillator potential. 

By the same Euler-Maclaurin method we have also investigated the anisotropic oscil- 
lator trap where the frequencies Lu x ,uj y and uj z are all different [[!(]] . Instead of the above 



8 



parameter b, we now introduce the three parameters = fiLOijk^T '. Corresponding to 
equation (|l3|), we then find for the first two leading terms in the result for the number of 
excited particles, 

n - - ^" 3(a> + + i + ss) Lhw + ■ ■ ■ <2i> 

Grossmann and HolthausQ considered a potential with lo x = 600 Hz, u> y = y2uj x and 
u z = \/3wx- They write the coefficient of Li2{\) as ^{ksT /%ui) 2 where uj = {uj x u y uj y ) 1 ^ 
and find 7 ~ 1.6 from a numerical summation over the energy levels of the anisotropic 
oscillator. On the other hand, we obtain from (|2l|) 



3X 1 / 3 / 1 1 1 



7= UJ U + V! + VSJ" L538 (22) 

which agrees quite well with their approximate result. 

We want to thank Mark Burgess for making referenceQ known to us. 
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